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ABSTRACT 



We update the list of GeV-TeV extragalactic y-ray sources using the 2-year catalog from the Fermi LAT and recent results ground- 
based y-ray telescopes. Breaks in the spectra between the high energy (100MeV< E <300GeV) and the very high energy 
(E> 200 GeV) ranges, and their dependence on distance, are discussed in the context of absorption on the extragalactic background 
light (EBL). We calculate the size of the expected break using a model for the EBL and compare it to the data taking into account 
systematic uncertainties in the measurements. We develop a novel Bayeasian model to describe this dataset and use it to constrain two 
simple models for the EBL-induced breaks. 



1. Introduction 

The extragalactic background light (EBL) is a diffuse field of 
U.V., optical and infra-red photons, with wavelengths in the 
range A = 0.1 -1000/mi, on which the integrated history of star 
formation in the Universe is imprinted. The spectral energy dis- 
tribution of the EBL consists of two distinct components: the 
first, peaking in vF v around ljum and commonly referred to as 
the cosmic optical background (COB), was produced by thermal 
emission from stars since the big bang. The second component, 
peaking at longer wavelengths (^ 100/^m), having comparable 
peak energy density to the COB and being referred to as the cos- 
mic infra-red background (CIB), origina tes from the absorption 
and reemission of starlight by dust (see lHauser & Dwekll200ll 
for review). 

Direct measurements of the EBL density are difficult due to 
local foregrounds, such as th e zodiacal light and Galactic radia- 
tion (lHauser & DweklfeOOll) . and are often interpreted as upper 
limits, while galaxy number counts in optical or infrared provide 
lower limits jMadau & Pozzettill2000l) . 

Since y rays of observed energy E y can interact with EBL 
photons of energy £ebl at a redshift z through yy — > e + + e~ 
when Ey/l TeV > 0.26 eV/£ , EBL(l + z), the spectra of distant 
extragalactic sources measured in the very high energy (VHE, 
E> 200 GeV) regime should differ from their emitted (intrin- 
sic) spectra if the EBL density is nonzercQ. Since a large frac- 
tion of the emitted power in BL Lac-type blazars is in y-ray 
band, this must be accounted for when spectral energy distri- 
butions are used to model their underlying physical properties 
dCoppi & Aharonianll999h . 

Finding clear evidence for this EBL-induced attenuation has 
proven remarkably difficult to date. The fall-off in the EBL spec- 
tral density between the COB and CIB peaks (around 0.1 eV) 
should be visible as a kink in the measured VHE spectra around 
1 TeV. This was sought for, e.g. in the blazar H 1426+428 by 



1 The cre ated pairs can a lso upscatter CMBR photons to high en- 
ergy y-rays (jprotheroe 1986) and induce yet an other spectral distortion 
mostly at energies ~ 100 GeV and below (e.g.. lAharonian etal 1 12002k 
iD'Avezac et alj2007l). a feature which has currently not jbeen observed 
(N eronov & Vovkl201(3h with instruments sensitive in that energy range. 



lAharonian et al. I d2003h . but results have been inconclusive given 
the large statistical errors. The signature of the EBL should also 
be evident in studies of the global population of VHE sources, 
since the energy-dependent attenuation increases with distance, 
such that the observed spectra are expected to become softer, i.e. 
the photon index, F, in power-law spectral fits should increase 
with redshift, z. Such studies have not been successful either, 
with no evidence for a redshift-depen dent effect being found by 
iMoril (l2009l) . lDe Angelis etail J2009L the authors attributing this 
to varying spect ral states inducing a large scatter in the data) or 
see in particular Figure 13). 
It has however been possible to constrain the EBL density 
in the energy range where they interact with observed y rays. 
Using VHE spectra from distant BL Lac objects and a theoret- 
ically motivated conjecture that the photon index of the intrin- 
sic spectrum cannot be harder than T/ a 1.5, several authors 
have derived upper li mits for the density close to the lower limits 
from galaxy co unts (|Dwek & Krennrichll2005t lAharonian et al 



2006). Recentlv. lAckermann et al.ld2012l) and lAbramowski et al 



d2013l) have measured the EBL density using its imprint in the 
spectra of BL Lac objects and found a density of EBL c ompati- 
ble with the best upper limits to date dMever et al.ll2.Q12b . 

Operating as a n all-sky monitor, the Fermi LAT 
dAtwood et alj 120091) observes y rays in the high energy 
(HE, lOOMeV < E < 300 GeV) range, where the effects of the 
EBL are much smaller than in the VHE. Sources detected both 
in the HE by Fermi and in the VHE by imaging atmospheric 
Cherenkov telescopes (IACTs) then provide an opportunity to 
probe the effects of the energy-dependent attenuation from the 
EBL absorption across a much wider energy range. Here we 
present an updated list o f Ge V- TeV sources, building on the 
work of lAbdo et all (120091) and lAckermann et al l (120111) . 



2. Selection of the sources 

Since the first detection of an extragalactic y-ray source in the 
VHE range by Whipple dPunch et all 19921) . 50 AGN have been 
discovered in this energy band. The rate of detections has in- 
creased dramatically with the increased sensitivity of the latest 
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generation of IACTs (VERITAS, MAGIC and H.E.S.S.) and an 
improved observation strategy using data from the Fermi -LAT. 
An up to date view of the VHE sky can be fo und by browsing 
the TeVCat catalogued (IWakelv & Horanll2008l) . 

We select AGN from TeVCat for which a HE and VHE spec- 
trum has been published and a firm redshift has been determined. 
Of the 50 AGN, 30 are BL Lacs with published spectral infor- 
mation. Three of these, namely 3C 66A, PKS 1424+240 and 
PG 1553+113, do not have a firm redshift determination. In addi- 
tion there are two VHE-detected FSRQs, 3C 279 and 4C +21 .35. 
Despite their large redshift (see Table [1), internal absorption 
close to the emission region may strongly affec t their spectra 
dAharonian et al. 2008; Sitarek & Bednarekl2008l) . hence we in- 
clude them in this study only for illustrative purposes. We also 
include the radio galaxies Centaurus A and M 87 in the sample. 

The second Fer mi catalogue of AGN (2LAC, 
lAckermann et al.l 1201 ll) includes 1057 sources associated 
with AGN of many kinds. In this list, 36 out of the 37 VHE 
BL Lacs have a Fermi counterpart, the six that are not detected 
being SHBL J00 1355.9-1 85406, 1ES 0229+200, 1ES 0347-121, 
PKS 0548-322, HESS J1943+213, 1ES 1312-423. 

By merging the two lists, our sample contains 23 HE- VHE 
BL Lacs and two radio galaxies. The characteristics of these 
sources are listed in Table Q] including the photon indexes from 
published power-law spectral fits in the HE (from 2LAC) and 
VHE (from reference given in the table), which form the dataset 
for this study. 

AGN are observed to be variable in all wavelengths. In the 
VHE regime, flaring episodes have been observed from a num- 
ber of such sources, in particular those that are bright and/or 
close by (such as Mrk 421, Mrk 501 and PKS 2155-304), but 
most have not shown clear evidence of variability - many have 
been detected close to the sensitivity threshold of the instrument. 
It is not improbable that variability is a common feature of HBLs 
in the VHE regime, future instruments will be able to probe this 
in a larger population of fainter, more distant sources. In the 
HE range, B L Lacs, and espec ially HBL, are found to be the 
less variable (lAbdo et al.ll2010ah . To reduce any bias introduced 
by the use of non-simultaneous observations, we use the VHE 
spectrum which has the lowest flux reported in the literature, re- 
sulting in a ge nerally good agre ement between the overlapping 
energy ranges (lAbdo et al.ll2009l 

3. Interpretation 

3. 1 . Spectral evolution with the redshift 

The mean HE and VHE indexes of our sample are (The) = 1 -86 
and (Tvhe) = 3.18, respectively. For each source the photon 
index measured in the VHE range is greater than or compatible 
with that found in the HE, i.e. AL = Fvhe - The ^ 0. In the HE 
band, the RMS of the measured indexes is (The = 0.26, and the 
excess variance, which accounts for the measurement errors^, is 
or™ E = 0.24. In the VHE regime, the RMS is o- VHE = 0.49, while 
the excess variance is o~y S HE = 0.10, showing that most of the 
sample variance can be ascribed to the errors on the individual 
measurements rather than to the intrinic distribution. 

The points on Figure \T\ show Ar versus the redshift z for our 
sample of sources. The two close-by radio galaxies do not show 
significant spectral breaks. For all other sources Ar > 0.5, and 

2 http://tevcat.uchicago.edu/ 

3 We define the excess variance of a set of measured quantities x, ± <r, 

as (cr^) 2 = {x 2 l )-{x,) 2 -(a 2 l ). 



those more distant than z — 0.1 exhibit a break of AF > 1.0. A 
dependence of Ar with z is apparent in our sample. 

If we were to assume that the intrinsic spectrum of each 
object was well represented by a single power law across the 
entire HE and VHE domain, as seems to be the case for the 
two nearby radio galaxies for which Ar ~ 0, we should expect 
that any significant break in the measured spectrum is the re- 
sult of absorption on the EBL, i.e. Fvhe = ^im + ATebl(E,z) ~ 
Ti„, + t/iogg C^' z )> where t(E 7 , z) is the optical depth due to the at- 
tenuation by pair production (lAbdo et al 1 I2009L Equation 2). To 
leading order the expected EBL break increases linearly in the 
redshift and we would AT to be correlated with z. The Pearson 
correlation factor for our dataset is p = 0.56 + 0.1 1, more than 
5cr away from 0, showing clear evidence that this dependency 
exists. 

In appendix IB. 21 we outline a simple method to evaluate 
the size of the break that might be expected from redshifting 
a curved intrinsic SSC spectrum within the HE and VHE ob- 
servation windows, which gives rise a K correction for the mea- 
sured spectral indexes. This study shows that such a K correction 
would account for ~ 15% of the observed break for the most dis- 
tant source in the sample at z — 0.3. 

3.2. Expected EBL-induced spectral break 

To further evaluate the data, we estimate the size of the spectral 
break as a function of re dshift from the EBL density model of 
iFranceschini et al.l (120081 hereafter Fra08). We perform a sim- 
ple simulation in which a hypothetical source with a flat spec- 
trum (Tj = 0) is placed at a distance z and its flux attenuated by 
the EBL. The simulated spectrum consists of 20 logarithmically 
spaced bins per decade, equally weighted, which is fitted with a 
power law model above 200 GeV to evaluate the measured in- 
dex. The limited photon flux at the highest energies is accounted 
for in an ad hoc manner; the upper energy bound of the fit is 
chosen to be the point at which the differential flux is 1% of the 
flux at 200 GeV (up to a maximum of 10 TeV). 

The predicted AF(z) is the black line shown in Figure Q] The 
shaded gray areas show uncertainties on this calculation, which 
arise from: 

- the ~ 10% energy resolution typical of IACTs which is taken 
into account by shifting the energy bins by +10% (dark gray 
area), 

- the threshold energy of the observations, which can vary 
from 100 GeV to more than 500 GeV (gray area), and 

- the systematic error on the measured photon spectral index, 
typically 0.2 (light gray). 

It is clear from Figure [1] that, for the majority of sources, 
the observed break, Ar, is systematically larger than that pre- 
dicted by the EBL model, ATebl- This is most notibly the case 
for 1ES 2344+514 (z = 0.044), PKS 2005-489 (z = 0.071), 
W Comae (z = 0.102), PKS 2155-304 (z = 0.116) and 
H 1426+428 (z = 0.129). The difference, AF - Ar £BL , is al- 
most certainly the result of (convex) intrinsic c urvature in the 
spect ra of these objects, which is not unexpected dPerlman et al.l 
120051) . and can have several interpretations, for instance as be- 
ing due to a turn-over in the distribution of the underlying em it- 
ting particles (acceleration effects; e.g.. iMassaro et al]l2006l) or 
to Klein-Nishina suppression (emission effects). 

Nevertheless, it is striking that there are many sources for 
which AF =i AT ebl- For these sources, the intrinsic broad-band 
y-ray spectra are compatible (within errors) with single power 
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laws. For those that additionally have The 5= 2.0, the high energy 
peaks are not constrained by the current observations, despite 
having a well defined observational vF v peak. The most striking 
examples are H 2356-309 (z = 0.129), 1RXS J101015.9-31 1909 
(z = 0.142), 1ES 1101-232 (z = 0.186), 1ES 0414+009 (z = 
0.287) and S5 0716+714 (z = 0.300). 

3.3. Constraining the EBL density 

Since a significant fraction of the observed break can be directly 
attributed to the EBL, we attempt to constrain its density by ap- 
plying a Bayesian model which takes into account the effects 
discussed in the previous section. The model is described fully 
in Appendix [A] We use two prescriptions to account for the ef- 
fects of the EBL on the spectra, ATebl{E,z). In the first, the 
break is mod eled as a linear function of the redshift, with coeffi- 
cient a, (as in lStecker & Scul ly 2006) and assume that the VHE 
measurements cover approximately the same energy range, so 
that the effect of energy threshold can be neglected. 

Ar £BL ,( fl ) = Ar EBL (Ei, Zi \a) * azi + O(z-). (1) 

In the second model we apply a scaling factor, a, to the EBL 
model of Fra08, which results in an expected break of, 

Ar £BL i(a) = ATebUEu Z/\a) = aAF Fra0& (E if z,), (2) 

where AFf ra o8(£;,ZiX is calculated for each source as in sec- 
tion [3?2] 

For both models, the posterior probability is computed and 
the results are given in Table [2] Figure |2]depicts the resulting Ar 
for each model, using the mean value (i.e. (a) and (a)) and the 
95% confidence level (CL) lower limit (i.e. ap<95% and ap<95%). 

The Bayesian model gives a value of (a) = 5.37 + 0.65 and 
an 95% CL upper limit of 6.44 for the linear EBL model, signif- 
icantly less than the value of 8.4+ 1.0 reported bv lYang & Wand 
(1201 Oh using a simple y 2 fit, which di d not account for the in- 
trinsic breaks. IStecker&ScuTi^(l2010l) found that their baseline 
model can be approximated by a linear coefficient of 7.99, which 
cannot be reconciled with the results presented here. The null hy- 
pothesis, i.e. that there is no dependence of AT with the redshift 
(a = 0), is rejected at more than 8cr. The spectral break predicted 
using the model of Fra08 is in good agreement with the data; the 
mean scaling factor is (a) = 0.85 + 0.10 and the 95% CL limit 
is a < 1.02. 

4. Conclusions and perspectives 

We have shown that broad-band y-ray spectra (from ~ 100 MeV 
to a few TeV) carry the imprint of the EBL and provide a unique 
dataset to probe its properties. The redshift dependence of the 
difference of the photon indices in VHE and HE range, AT, 
is found to be compatible with expectations from EBL atten- 
uation. The Pearson correlation coefficient shows that Ar and 
z are significantly correlated. We developed a Bayesian model 
to fit the data set, accounting for intrinsic spectral softening, 
and find that the EBL densi t y is co nsistent with the value pre- 
di cted by iFranceschini e t al. (2008 1). Similar results were found 
bv lAckermann et al.l (120121) and I Abramowski et al.l (120131) . who 
modeled the EBL-absorbed spectra of AGNs detected in the 
HE and VHE regimes respectively, and found scaling factors of 
o-Fermi = 1 -02 + 0.23 and (Xhess = 1-27*5 J?. Their approach has 
the potential to be more powerful than that used here, since it can 
probe the features of the EBL-absorption signature as a function 
of energy. However their approach is more reliant on the detailed 



modeling of the intrinsic spectra of the objects and of the EBL 
density and does not take advantage of the wider energy band 
available when the HE and VHE observations are combined. The 
two approaches are complimentary and yield roughly compati- 
ble results within the combined statistical and systematic errors. 
Taking only the statistical errors, a^ 2 fit to the HESS, Fermi and 
our results gives a mean value of (X C ombined ~ 0.98 and a value 
of x 1 = 5.45 for two degrees of freedom, compatible with the 
hypothesis that the values are consistent at the 1.85cr level. Our 
model also offers a simple prescription for constraining the red- 
shift of a GeV-TeV sources based on their measured value of AT 
(see Appendix lAl. Applying our findings to PG 1553+113 and 
3C 66A leads to z < 0.64 and z < 0.55 respectively, in good 
agreement with the spectroscopic constraints. 

No sources have breaks significantly smaller than AFebl- 
Significant deviations from the expected EBL-induced spectral 
breaks could indicate either concave curvature in the intrinsic 
spectrum of the source or that other processes are at play during 
the propagation of y rays, such as cosmic-ray interactions along 
the line of sight which create spectral so ftening in high-redshift 
source spectra (lEssev & Kusenko1l2012l) . Our findings indicate 
that experimental uncertainties need to improve, and firmer red- 
shift estimations established, before the significance of this ef- 
fect can be assessed. 

The recent commissioning of the 28m H.E.S.S. 2 telescope, 
the commissioning of the upgraded MAGIC telescopes and the 
upgrade of the VERITAS cameras and trigger, should increase 
the distance at whic h new blazars ca n be detected, while the 
planned CTA project (lActis et al.l201 ll) should reduce the uncer- 
tainties mentioned above due to its superior sensitivity. Finally, 
a useful feature of studies of spectral breaks vs redshift, such as 
this one, is their capacity to provide distance estimations for BL 
Lacs (see, e.g., lAbdo et al.lfeoiObl: iPrandini et aDl2010i 1201 2l) 
since an estimated 50% of this population has an unknown or 
uncertain redshift. 



References 

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010a, ApJ, 722, 520 
Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, ApJ, 707, 1310 
Abdo, A. A., Ackermami, M., Ajello, M., et al. 2010b, ApJ, 708, 1310 
Abramowski, A., Acero, F., Aharonian, F, et al. 2013, A&A, 550, A4 
Abramowski, A., Acero, F, Aharonian, F, et al. 2012, A&A, 542, A94 
Abramowski, A., Acero, F, Aharonian, F., et al. 2010, A&A, 516, A56 
Acciari, V. A., Aliu, E., Arlen, T., et al. 2010, ApJ, 715, L49 
Acero, F, Aharonian, F, Akhperjanian, A. G., et al. 2010, A&A, 511, A52 
Ackermann, M., Ajello, M., Allafort, A., et al. 2012, Science, 338, 1190 
Ackermann, M. et al. 201 1, ApJ, 743, 171 

Actis, M., Agnetta, G., Aharonian. F., et al. 201 1, Experimental Astronomy, 32, 
193 

Aharonian, F, Akhperjanian, A., Beilicke, M., et al. 2003, A&A, 403, 523 
Aharonian, F, Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, Nature, 440, 
1018 

Aharonian, F. A., Khangulyan, D., & Costamante. L. 2008, MNRAS, 387, 1206 
Aharonian, F. A., Timokhin, A. N., & Plyasheshnikov, A. V. 2002, A&A, 384, 
834 

Ahrendt, P. 2005, The Multivariate Gaussian Probability Distribution, Tech. Rep. 
IMM3312, DTU 

Aleksic, J., Alvarez, E. A., Antonelli, L. A., et al. 2012a, A&A, 539, A118 

Aleksic, J., Alvarez, E. A., Antonelli, L. A., et al. 2012b, ApJ, 748, 46 

Aleksic, J., Antonelli, L. A., Antoranz, P., et al. 201 1, ApJ, 730, L8 

Aliu, E., Aune, T, Beilicke, M., et al. 201 1, ApJ, 742, 127 

Anderhub, H., Antonelli, L. A., Antoranz, P., et al. 2009, ApJ, 704, L129 

Atwood, W. B. et al. 2009, ApJ, 697, 1071 

Band, D. L. & Grindlay, J. E. 1985, ApJ, 298, 128 

Benbow, W. & The VERITAS Collaboration. 2011, in Proc 32nd Intl. Cosmic 

Ray Conf., Beijing, arXiv e-prints, 1 1 10.0038 
Coppi, P. S. & Aharonian, F. A. 1999, Astroparticle Physics, 1 1, 35 
D'Avezac, P., Dubus, G., & Giebels, B. 2007, A&A, 469, 857 



3 



D.A. Sanchez et al.: Evidence for a cosmological effect in y-ray spectra of BL Lacs. 




j ! i i i i ! i i i i ! i i i i ! i i i i ! i i i i ! i i i i ! I i_ 

0.1 0.2 0.3 0.4 0.5 0.6 

Redshift z 

Fig. 1. The value of AF as a function of the redshift z- The black line is the theoretical break obtained with iFranceschini et al.1 
(2008) model. The uncertainties due to the energy resolution of IACTs (dark gray), the different threshold energies (gray) and the 
systematic errors of 0.2 (light gray) are shown. The FSRQs and the BL Lacs with uncertain redshift are shown for illustration. 
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Fig. 2. The value of AF as a function of the redshift z. Overlaid are the predicted breaks obtained from the Bayesian fit (simple line), 
as well as the 95% CL lower limits (hashed area). Left is the linear approximation (Equation[T) and right is the scaled model of 
Fra08 (Equation For both models, none of the sources are significantly lower than the 95% CL lower limit. The symbols are 
descibed in FigureQ] 
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Appendix A: Bayesian model 

A.I. Development of the model 

To extract information about the EBL density from the data 
presented in Figure Q] we devel op a "hierarchical" Bayesian 
model (using the terminology of iGelman et al.ll2003l) which is 
described by source-by-source spectral parameters and global 
parameters specifying, amongst other things, the EBL density of 
primary interest here. We use a Bayesian methodology to write 
the posterior density for the model parameters, marginalize over 
the source-by-source parameters which are not of interest and 
produce estimates and confidence intervals for the EBL density. 

In the development below we make repeated use of the con- 
ditional probability rule (CPR), that a joint probability of two 
(sets of) events A and B, P(A, B), can be expressed as P(A, B) = 
P(A\B)P(B). In the case that the two events are independent 
this becomes P(A, B) = P(A)P(B). We also frequently use the 
rule for marginalizing (or integrating) over unwanted parame- 
ters, P(A) = J P(A, B)dB. Finally, we us e the standard identity 
for the product of two Gaussians (see e.g. lAhrendtl2005b . In par- 
ticular if N(x\jj, a 2 ) denotes a Gaussian distribution with mean p 
and variance <x 2 , then, 



f 



N(x\/xu(rl)N(x\n2, o\)dx 



^{o- 2 l+ o- 2 2 y 



exp < - 



2{a\ + o- 2 2 ) 



(A.l) 



The data set to be modeled consists of N GeV and Te V spec- 
tral measurements, and with measurement variances of 
cr?,. and cn^. and redshifts z, (which are themselves not consid- 
ered as measurement data). In what follows we refer frequently 
to the measured spectral break, AFf = F^ - T^?.. We write the 
data set as Y = {Tg, 1^}. 



Each source is parameterized by four values, the intrinsic 
spectral indexes in the GeV and TeV regimes, F' Gj and FL and 
the spectral indexes after absorption by the EBL, rd. and Fj.. In 
what follows we refer frequently to the intrinsic spectral break, 
Ar' = F' Tj - FL.. The global parameters that describe the EBL 
absorption itself are denoted abstractly as G and we write the set 
of all parameters of the model as = {F Gi , F^., F Gl , Fj., G}. 

Bayes' theorem allows us to write the posterior probability 
distribution of the parameters after the measurements have been 
made, P(®\Y), in terms of the standard Likelihood of the data, 
P(Y\&), and the prior probability distribution of the model pa- 
rameters, P(0): 

P(&\Y) cc P(&)P(Y\&) 

The relation is written as a proportionality, instead of as an 
equality, as a global normalization factor has been neglected. 

The model has four primary components that we discuss be- 
low. These are: (1) the likelihood for the GeV and TeV measure- 
ments, given the true absorbed spectral indexes of the sources, 
(2) a relationship between the absorbed index in the GeV regime 
and the intrinsic (unabsorbed) index, (3) an relationship for the 
analogous indexes in the TeV regime, and (4) a specification of 
how the intrinsic index in the GeV regime is related to the intrin- 
sic index in the TeV regime. 

We assume the measurements of the individual indexes are 
independent (no correlation between measurements from differ- 
ent sources or between the GeV and TeV bands) and that the 
distribution for each measurement (F Gi or F^) is Gaussian with 
mean given by the appropriate absorbed index and with variance 
given by the measurement errors squared. Therefore the likeli- 
hood is, 

p(y\@) = n^l r G/^G,-W(r^ir^,,^.). 

i 

Similarly for the prior, we assume that the only link between 
the source-by-source index parameters for different sources 
comes through the global parameters G. Therefore, using the 
CPR we can write, 

p(&) = P(G) Y] P(r 7 G ,., r 7 ,., i* , f^ig) 

i 

It now remains only to describe how the intrinsic and absorbed 
indexes for each source are related. We assume that the absorbed 
index in each band depends only on the intrinsic index in that 
band and on the EBL parameters. Repeatedly applying the CPR, 
this gives, 

PPou Tt,> 4, 4|G) = P(T A Gi \r' Gi , G)P(F A i \F I T „ G)P(F ! T „ F 7 G ,.|G). 

(A.2) 

We further assume that there is no absorption in the GeV regime, 
and that the absorption in the TeV regime changes the index in 
a deterministic way. Specifically we assume that, (i) F Gi = F Gi 
and (ii) F A . = F' Ti + AF EBL (Ei,Zi\G), where AF EBL (Ej,Zi\G) is a 
function giving the change in TeV index for a source at redshift 
Zi measured at a TeV "threshold" energy of E,. This can be ex- 
pressed in terms of a probability using the Dirac (5-function: 

P(Xt i \r 1 Gi ,G) = s(F A ii -F Gi ) ® 

PiF^F'^, G) = S(F A - F' Ti - AF EBL (E h Zi \G)) (ii) 



AFEBLiEi^AGVPiF^FcP). 



Equation lA.2l then reads: 

PCr^.r^.r^.r^iG) 

= 6(F A Gi -F' Gi )6(F A Ti -F ! Ti 
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The final and most interesting part of the model is to describe 
how the intrinsic indexes in the two bands are related. We expand 
this using the CPR to give, 

p(r^,r 7 G ,.|G) = P(r^|r^,G)P(r^|G). 

We adopt a uniform prior for r^., P(TL.) = 1, and restict our- 
selves to forms for the conditional probability that can be ex- 
pressed as a function of the intrinsic break, P(r' T JTL., G) = 
P(AF[|G), and assume that the break for each source is drawn 
from a single universal distribution which has, at most, some de- 
pendence on the global parameter set, G. The final expression 
for the prior for the parameters of each source is: 



^(r^,r£,,r*|G) 

= 5(^-1^(^-1^ 



Ar £2>L (£ />Z/ |G))P(Ar[|G). 



Putting everything together, and recognizing that we haven't 
discussed the parameters G yet, and leaving the exact choice of 
prior for the intrinsic break open for the present time, the full 
posterior density is, 

^{r^r^,r^,r^.},G|{r«,r^})<x 

p(G) Y] P(MtiG)iV(rg.|r£ il c4)iv(r$n*, a 2 ^ 

i 

sKi - r'c^K - r x/ - Ar EBL (E hZi \G)) 

We marginalize over the parameters that are not of interest, 
r' Gi , r^ ; , rd. and to give the posterior probability for the 
global parameters G. The integrals over the parameters for each 
source can be done separately, 

h= J dT' a j dT' Ti j dT A Ci j dT* 

P(AT[|G)AT(r^.|r^, ar^NiT^, o*,) 

8<Tai ~ - Tt/ - ^EBL(E t ,Zi\G)) 

The integrals over rd. and r£. can be done immediately against 
the delta functions to give, 



J, 



, = j dT^j dT^PiAT^G) 

WaFoi' ^G/W(r^|r 7 ,. + AT EBL (E,, Zi \G), cr^) 



Making a change of integration variable from r^, ; to ATj, the 
Gaussians can be manipulated to give, 

I, = j d(AT[)P(AT[|G) J dr' a 

Wo^Gi'^WGiK - AT 7 - Ar EBL (E hZi \G),^. 



The second integral can be evaluated using equation lA.il Putting 
all the source integrals together gives the general expression, 



p(G|{r^,r^})<xP(G) 
Y] J d(Ar[)P(Ar[|G)Ar(Ar[|Arf-Ar £BL (E i)Zi |G),4 i +^ i ) 



A.2. Various priors for AT' 

We examine three concrete cases for P(Ar'|G). The first two are 
based on simple assumptions and result in analytic expressions 
for the full posterior probability that are easy to understand. The 
final prior distribution is derived from Monte Carlo realizations 
of an SSC model, as described in appendix lB.il The results from 
this final case that are presented in section l3~3l Here we develop 
the first two cases. A relatively weak assumption is that the TeV 
index can be no harder than the GeV index. It would seem rea- 
sonable to express this using the Heaviside step function, ©(x), 



P(ATflG) = ©(Aft). 



(A.4) 



However this is unsatisfactory, as it asserts that the mean intrin- 
sic break is infinite, (AT') — > oo, which is clearly not realistic. 
As will be seen below, this results in an unphysical posterior 
distribution. To remedy this failing we instead assume that the 
prior distribution of the intrinsic break is given by a Gaussian, 
truncated at negative values: 



p(Arj|G) = ®(Ar!)N(Ar!\m,oj). 



(A.5) 



In this case the prior is parameterized by two values, yU/ and 07, 
which must be either estimated in the problem (i.e. added to the 
global parameter set G), or specified externally. We will simply 
assume fij — and derive the results for various reasonable val- 
ues of 07. 

Equation lA.3l can be used to calculate the posterior distribu- 
tion in the three cases for P(AF') discussed above. In the second 
case, with the prior given by Equation IA. 51 the final integration 
can be done to give, after applying the formula for the product 
of Gaussians, 



p(G|{r£.,rE}, A i / = o,<r / )oc 



P(G)\ 



exp 



•M\2 



{AY EBL {E u Zi\G) - ATf) 
2(cr 2 Gi + <r\ + oj) 



erfc 



( ( r / (Ar £BL (£,,z,|G)-Arf) ) 

S i [ r t > ( A - 6 ) 

(V^o* +<7*)J(< + 4 i + °/) I J 



where erfc(;t:) = -4= J x dx' e is the complementary error 
function. 

It is instructive to examine the two limiting cases of 07 = 
and 07 — > 00. In the first case we arrive at, 



logP(GRT« r$},ju 7 = 0,cr 7 = 0) 



Z 



M *2 



(Ar £flL (£,-,z,|G)- Arf) 



which is exactly the expression that would result from a simple 
least-squares fit to the measured spectral breaks. In the second 
case (07 — > 00) we have, 



p(Gnrg.,r^},/i/ 



0, 07 — > 00) 



"] erfc j 



{ATEBLiEuZiW-ATl 



M 



(A.7) 



which is the same expression as would be derived starting from 
the prior given by Equation lA.4l We therefore have an expression 
that transforms continuously between the two clearly identifiable 
(A. 3) extremities as a function of 07. 
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As described in section [3731 we attempt to derive constraints 
on two simple models for the EBL. First, a linear approximation, 
given by Equation]]"] and second, a scaling of the EBL model of 
iFranceschini et al] ( 2008h . as described by Equation [2] In both 
cases we assume a uniform prior in the positive region of the 
parameter space for a, P(a) = 0(a), and a, P(a) = ®(a), respec- 
tively. 

Using either of these, the deficiencies in the model given by 
Equation lA.4l is finally evident. Combining Equations l2l and lA~7l 

we get, 

P(a\{T^,T^ i \, Ml = 0,o- 1 ^oo)cK 

0(a) erfc<^ — — — \. (A.8) 

The most probable value occurs at a = 0, which is consistent 
with the assertion in this case that (Ar') — > oo. 

The results in section 13.31 are derived from combining 
Equations Q] (or |2]i and IA.3I with the prior calculated in ap- 
pendix iBTl and integrating numerically. 



A.3. Results with different priors 



Table IA.1I presents the results obtained with the half-Gaussian 
prior of Equation IA.5I (with /i/ = 0), for four values for 
the variance of the distribution of the intrinsic break, 07 = 
0.25, 0.5, 1.0 and 2.0, and using the prior derived from SSC mod- 
eling, illustrated in Figure IbTTI 

With increasing value of 07, the most probable value of a or 
a decreases, since a lower EBL-induced break is needed to re- 
produce the data. Results derived with values of 07 = 0.5 and 1 .0 
are in good agreement with those from the SSC model. 

A.4. Constraints on the redshift 

The Bayesian methodology can also be used to constrain the red- 
shifts of GeV-TeV blazars from their measured spectral breaks. 
In the case of a single source, equations lA.3l andl2lcan be adopted 
to express the posterior probability for the parameters G = {a, z], 
given their priors P(a) and P(zfl This can then be marginalized 
over a to give, 

P(z\T%,r»!)ocP(z) j daP(a) j £ /(AF / )P(AF / )x 

jV(Ar'|Ar M - aAT Fnm (E, z\G), cr| + cr 2 ). (A.9) 

The prior for the redshift, P(z), could be estimated from the red- 
shift distribution of detected GeV-TeV detected blazars, i.e. from 
the values presented in Table Q] However the true distribution is 
probably not well represented by the small number of sources 
detected, so we seek an alternative approach. Another option is 
to use a flat distribution, which is conservative but also unreal- 
istic. We instead compromise and use the distribution of 2FGL 
BL Lac objects and unknown AGN. Since Fermi AGN are de- 
tected to higher redshifts than those at TeV energies, this should 
still be a conservative approach. For the prior on the EBL scal- 
ing, P(a), we use the positive portion of a Gaussian with mean 
1.0 and RMS 0.3, P(a) = 0(a)A^(a|l.O,O.3 2 ). Finally, we use 
the prior on Ar' derived from our SSC simulations, as described 
above. 



4 We neglect the dependence of the prior for z on the strength of the 
EBL (a), i.e. we assume incorrectly that P(z\a) = P(z). 



The mean redshift and upper limits for some values of AF, 
calculated using a typical value of <tL + a\ - 0.2 2 , are given in 
Table IA.2I The corresponding relation between (z) and AT can 
be approximated by 

<Z> ~ 0.024 + 0.079Ar + 0.022Ar 2 - O.OOlOAr 3 

and the relation between Zp<gs% and Ar by 

z P<95% * 0.081 +0.081Ar + 0.080Ar 2 - 0.01 1AF 3 . 

Appendix B: Synchrotron Self-Compton 
Simulations 

B.1. Determination of the intrinsic break properties 

In order to derive a plausible prior probability density for the 
intrinsic break between HE and VHE, for use in the Bayesian 
model, we produce a set of Monte Carlo simulations of hypo- 
thetical BL L acs using a one-zone sy nchrotron self-Compton 
(SSC) model dBand & Grindlavlll985l) . which is often used to 
successfully reproduce the time-average SED from radio to TeV 
energies. 

We assume a spherical emission zone, with a size R, moving 
at a bulk Doppler factor 6. This region is filled by a uniform mag- 
netic field B, and a population of electrons with a density N e (~y) is 
responsible for the synchrotron emission. The synchrotron pho- 
tons are upscattered by the same population of electrons to pro- 
duce y rays. 

The distribution of elect rons is described by a power-law 
with an exponential cut-off dLefa et al.ll2012l) . N e (y) = N y p ■ 
exp(-y/y cut ). The model therefore has three parameters to de- 
scribe the electron population (No, p and y cut ) and four to de- 
scribe the jet properties (z, R, 5 and B). Among these param- 
eters, R, 6 and No have only an achromatic effect, and z pro- 
duces a small K correction, which is evaluated separately in ap- 
pendix lB~2l The spectral break is also insensitive to the value of 
the index p of the electron distribution. The parameters which 
determine the intrinsic break are therefore B and y cut . 

We perform 10 5 simulations in which the values of B and 
y cut and uniformly drawn in the range 0.01 < B < 0.5 G and 
3 x 10 4 < y cu t < 1 x 10 7 . The other parameters are kept fixed at 
the values given in Table IB. fl which are typical for BL Lacs. 

Since only the BL Lac-type SEDs are of interest in this study, 
simulations for which the synchrotron emission peaks at ener- 
gies lower than 10 eV or higher that 10 MeV and having a HE 
index of F > 2 have been removed. The distribution of the in- 
trinsic break Fy HE - r^ E , depicted in Figure IbTTI has a sharp rise 
below 0.2 and a long tail at higher break values. 

B.2. Effects of the energy shift due to the distance 

SSC simulations can be used to evaluate the size of the K cor- 
rection required to account for the redshifting of the intrinsic 
spectrum into the fixed HE and VHE observation windows. This 
effect produces a trend of increasing observed AT with z, even in 
the absence of EBL absorption. 

As before, the SSC parameters are fixed to the values in 
Table |B7T] but with B = 0.1 G and y cut = 1.6 X 10 5 . This pro- 
duces an intrinsic spectrum with r^ E = 1.85 and Ty HE = 3.04, 
i.e. AFj = 1.19. Simulating the same source with redshifts in 
the range 10~ 4 < z < 0.7, and with no EBL absorption, leads 
to an additional component in the observed Ar which increases 
with redshift, as depicted in Figure lB~2l This corresponds to the 
K correction for this intrinsic spectral shape and is too small to 
explain the trend observed in the data. 
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Table B.l. Parameters used for the SSC simulations. 



Parameters 


Value 


z 


0.1 


R 


4.5 x 10 16 cm 


6 


20 


P 


2.3 


No 


3 x 10 3 cm" 3 


B 


0.01- 0.5 G 


7cut 


3 x 10 4 - 1 x 10 7 
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Table 1. List of HE-VHE BL Lacs and radio galaxies. Only statistical errors are given. 



Source Name 


Q'KOOO 


<5j2000 


Type 


- 




r 


/HE 


Ref. 


Centaums A 


13 h 25 m 27.6 s 


-43° 01' 09" 


FR1 


0.00183 


2.76 ± 


0.05 


2.7 


±0.5 


[1] 


M 87 


12 h 30 m 49.4 s 


+ 12° 23' 28" 


FR1 


0.004233 


2.17 ± 


0.07 


2.60 


+ 0.35 


[1] 


Markarian 421 


ll h 04 m 27.3 s 


+38° 12' 32" 


HBL 


0.031 


1.77 ± 


0.01 


2.20 


±0.08 


[1] 


Markarian 501 


16 h 53 m 52.2 s 


+39° 45' 37" 


HBL 


0.034 


1.74 ± 


0.03 


2.54 


±0.70 


[1] 


1ES 2344+514 


23 h 47 m 04.8 s 


+51° 42' 18" 


HBL 


0.044 


1.72 + 


0.08 


2.95 


±0.12 


[1] 


Markarian 180 


ll h 36 m 26.4 s 


+70° 09' 27" 


HBL 


0.046 


1.74 ± 


0.08 


3.3 


±0.7 


[1] 


1ES 1959+650 


19 h 59 m 59.9 s 


+65° 08' 55" 


HBL 


0.048 


1.94 + 


0.03 


2.58 


±0.18 


[1] 


BL Lacertae 


22 h 02 m 43.3 s 


+42° 16' 40" 


LBL 


0.069 


2.11 ± 


0.04 


3.6 


±0.6 


[1] 


PKS 2005-489 


20 h 09 m 25.4 s 


-48° 49' 54" 


HBL 


0.071 


1.78 ± 


0.05 


3.20 


±0.16 


[8] 


RGB J0152+017 


01 h 52 m 39.6 s 


+01° 47' 17" 


HBL 


0.080 


1.79 + 


0.14 


2.95 


±0.36 


[1] 


W Comae 


12 h 21 m 31.7 s 


+28° 13' 59" 


IBL 


0.102 


2.02 + 


0.03 


3.81 


±0.35 


[1] 


PKS 2155-304 


21 h 58 m 52.1 s 


-30° 13' 32" 


HBL 


0.117 


1.84 ± 


0.02 


3.32 


±0.06 


[1] 


B3 2247+381 


22 h 50 m 06.6 s 


+38° 25' 58" 


HBL 


0.119 


1.83 + 


0.11 


3.2 


±0.6 


[9] 


RGB J0710+591 


07 h 10 m 30.1 s 


+59° 08' 20" 


HBL 


0.125 


1.53 ± 


0.12 


2.69 


±0.26 


[5] 


H 1426+428 


14 h 28 m 32.7 s 


+42° 40' 21" 


HBL 


0.129 


1.32 ± 


0.12 


3.5: 


t0.35 


[1] 


1ES 0806+524 


08 h 09 m 49.2 s 


+52° 18' 58" 


HBL 


0.138 


1.94 ± 


0.06 


3.6 


± 1.0 


[1] 


1RXS J101015.9-311909 


10 11 10 m 15.03 s 


-31° 18' 18.4" 


HBL 


0.142 


2.09 ± 


0.15 


3.08 


± .42 


[7] 


H 2356-309 


23 h 59 m 07.9 s 


-30° 37' 41" 


HBL 


0.167 


1.89 + 


0.17 


3.06 


± 0.15 


[10] 


RX J0648.7+1516 


06 h 48 m 45.6 s 


+ 15° 16' 12" 


HBL 


0.179 


1.74 ± 


0.11 


4.4 


±0.8 


[4] 


1ES 1218+304 


12 h 21 m 21.9 s 


+30° 10' 37" 


HBL 


0.182 


1.71 ± 


0.07 


3.08 


±0.34 


[1] 


1ES 1101-232 


ll h 03 m 37.6 s 


-23° 29' 30" 


HBL 


0.186 


1.80 + 


0.21 


2.94 


±0.20 


[1] 


RBS 0413 


03 h 19 m 51.8 s 


+ 18° 45' 34" 


HBL 


0.19 


1.55 ± 


0.11 


3.18 


±0.68 


[2] 


1ES 1011+496 


10 h 15 m 04.1 s 


+49° 26' 01" 


HBL 


0.212 


1.72 + 


0.04 


4.0 


±0.5 


[1] 


1ES 0414+009 


04 h 16 m 52.4 s 


+01° 05' 24" 


HBL 


0.287 


1.98 ± 


0.16 


3.45 


±0.25 


[3] 


S5 0716+714 


07 h 21 m 53.4 s 


+71° 20' 36" 


LBL 


0.300 


2.00 ± 


0.02 


3.45 


±0.54 


[6] 



Additional sources used for illustration only: 

3C66A 02 h 22 m 39.6 s +43° 02' 08" IBL 0444? 1.85 ±0.02 4.1 ±0.4 [IT 
4C +21.35 12 h 24 m 54.4 s +21° 22' 46" FSRQ 0.432 1.95 ±0.21 3.75 ± 0.27 [11] 
PG 1553+113 15 h 55 m 43.0 s +11°11'24" HBL 0.43 -0.58 1.67±0.02 4.41 ±0.14 f [12] 
3C 279 12 h 56 m 11.2 s -05° 47' 22" FSRQ 0.536 2.22 ± 0.02 4.1 ±0.7 [1] 



Re ferences. [11 lAbdo et al.1 d2009l. see refe renc es therein); [2] IVERITAS Collaboration et all d2012l) ; f31 IH.E.S.S. Collaboration et alj d2012h; 
[4llAliuetal.U201lh: r5llAcciari et al.ld20irl: r61 r Anderhub et al.U2009]): mlAbramowski et alJd2012h: mlAcero et al.ld2010h : r9l lAleksic et al.1 
d2012ah; riOHAbramowski et al.1 d2010h ; llll lAleksic et al.ld2011h : ri21 lBenbow & The VERITAS Collaboration! d201 ll) . 
t lAleksic et al.1 d2012bh recently published long-term observations on this object and derived a compatible photon index. 

Table 2. Summary of results with the linear parametrization and scaled iFranceschini et al.l d2008l) EBL model obtained with the 
Bayesian approach described in the text. 



Parameter 


Linear model 


Parameter 


Scaled Fra08 model 


Mean value: (a) 


5.37 


Mean value: (a) 


0.85 


RMS: V<« 2 > " (a) 1 


0.65 


RMS: V<a 2 >-<<*> 2 


0.10 


Upper limit: a P< g$% 


6.44 


Upper limit: a P< 95% 


1.02 


Lower limit: a P>S c /o 


4.32 


Lower limit: a P>5 % 


0.69 



Table A.l. Summary of results from Bayesian model with different priors. 



Parameter 


Half-Gaussian prior (EquationlA.51 /i/ = 0) 
a, = 0.25 cr, = 0.5 a, = 1 a, = 2.0 


SSC prior 
(Figure ED 


Mean value: (a) 


7.40 


6.04 


4.64 


3.12 


5.37 


RMS: yj(a 2 )-(a) 2 


0.55 


0.64 


0.84 


1.24 


0.65 


Upper limit: a P< g=,% 


8.29 


7.08 


5.97 


5.02 


6.4 


Lower limit: a P>s % 


6.5 


4.99 


3.20 


0.9 


4.32 


Mean value: (a) 


1.16 


0.96 


0.74 


0.50 


0.85 


RMS: yj(a 2 )-(a) 2 


0.08 


0.10 


0.13 


0.20 


0.10 


Upper limit: a P< ^% 


1.31 


1.12 


0.95 


0.80 


1.02 


Lower limit: a P>5 % 


1.03 


0.80 


0.52 


0.15 


0.69 
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Table A.2. Prediction of mean redshift value (z) and upper limit zp<95% for a give value of AT. 



AT 


<z> 


Zp<95% 


0.5 


0.07 


0.14 


1.0 


0.12 


0.24 


1.5 


0.19 


0.34 


2.0 


0.26 


0.46 


2.5 


0.34 


0.60 


3.0 


0.43 


0.75 


3.5 


0.53 


0.86 


4.0 


0.62 


0.95 



< 

CL 
T3 



0.06 



W 0.05 

a> 
> 

= 0.04 

!5 

-Q 

2 0.03 



0.02 



0.01 



~i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — r 




J U_l I I I I I I I I I I I I I I I I I I I I I I I I I I I L 



0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 

^TeV'^GeV 

Fig. B.l. Tvhe - The distribution (red histogram) obtained with the SSC simulations and application of cuts described in the text. 
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